#include "KSADAnalyze.hh"
#include "KSADAnalyzeVis.hh"
#include "KSADAnalyzeMessenger.hh"
#include "KSADAnalyzeEvent.hh"
#include "KSADPrimaryGeneratorAction.hh"
#include "G4RunManager.hh"
#include "KSADGlobals.hh"
#include "TFile.h"
#include "TApplication.h"
#include "TH3.h"
#include "TH2.h"
#include <iostream>
#include <fstream>

KSADAnalyze::KSADAnalyze()
{
	eventList = new KSADAnalyzeEvent*[0];
	numberOfEvents = 0;
	messenger = new KSADAnalyzeMessenger(this);
	h3 = 0;
	maxValue = 0;
}

KSADAnalyze::~KSADAnalyze()
{
	for(int i = 0;i < numberOfEvents;i++) delete eventList[i];
	delete eventList;

	if(h3) h3->Delete();
	if(gFile) gFile->Close();
}
	
void KSADAnalyze::Analyze()
{
	if(dataFilePath != "")
	{
		if(gFile) gFile->Close();
		gFile->Open(dataFilePath);
		h3 = (TH3D*)(gFile->Get("KSADh3"));
	}
	else h3 = globalH3;

	if(h3)
	{
		maxValue = h3->GetMaximum()*MeV;
		numberOfEvents = IdentifyEvents();
	}
	else G4cout << "No h3\n";

	G4cout << numberOfEvents << G4endl;
}

int KSADAnalyze::IdentifyEvents()
{
	if(h3->GetBinContent(h3->GetNbinsX()+1,0,0) > 0)
	{
		G4cout << "Warning: data histogram has hits in overflow, aborting event identification.\n";
		return 0;
	}

	int numEvents = 0;
	int **vetoList = new int*[0];
	int vetoCount = 0;

	for(int l = 0;l < numberOfEvents;l++) delete eventList[l];
	delete eventList;

	eventList = new KSADAnalyzeEvent*[0];

	for(int i = 0;i < h3->GetNbinsX();i++)
	{
		for(int j = 0;j < h3->GetNbinsY();j++)
		{
			for(int t = 0;t < h3->GetNbinsZ();t++)
			{
				// If condition is satisfied, event is found, and the process to add
				// the current crystal's surrounding energy deposits is started
				if(!InList(vetoList,vetoCount,i,j,t) && h3->GetBinContent(i+1,j+1,t+1) > 0.000100)
				{
					numEvents++;
					int iLow = i; int jLow = j; int iHigh = i; int jHigh = j; int tLow = t; int tHigh = t;
					vetoList = AddToList(vetoList,vetoCount,i,j,t);
	
					KSADAnalyzeEvent** eventListOld = eventList;
					eventList = new KSADAnalyzeEvent*[numEvents];

					for(int l = 0;l < numEvents-1;l++)
					{
						eventList[l] = eventListOld[l];
					}

					delete eventListOld;
					eventList[numEvents-1] = new KSADAnalyzeEvent();
	
					int** groupVetoList = new int*[0];
					int groupVetoCount = 0;
					groupVetoList = AddToList(groupVetoList,groupVetoCount,i,j,t);
	
					int newVetos = 1;

					// Find surrounding energy deposits
					while(newVetos > 0)
					{
						int iNew = groupVetoList[groupVetoCount-newVetos][0];
						int jNew = groupVetoList[groupVetoCount-newVetos][1];
						int tNew = groupVetoList[groupVetoCount-newVetos][2];
	
						for(int di = -2;di < 3;di++)
						{
							for(int dj = -2;dj < 3;dj++)
							{
								for(int dt = -2;dt < 3;dt++)
								{
									int feelerI = iNew+di;
									int feelerJ = jNew+dj;
									int feelerT = tNew+dt;
									double feelerContent = h3->GetBinContent(feelerI+1,feelerJ+1,feelerT+1);

									if(feelerContent > 0)
									{
										if(!InList(vetoList,vetoCount,iNew+di,jNew+dj,tNew+dt))
										{
											vetoList = AddToList(vetoList,vetoCount,iNew+di,jNew+dj,tNew+dt);
											groupVetoList = AddToList(groupVetoList,groupVetoCount,iNew+di,jNew+dj,tNew+dt);
											newVetos++;

											if(iLow > iNew) iLow = iNew;
											else if(iHigh < iNew) iHigh = iNew;
	
											if(jLow > jNew) jLow = jNew;
											else if(jHigh < jNew) jHigh = jNew;
	
											if(tLow > tNew) tLow = tNew;
											else if(tHigh < tNew) tHigh = tNew;
										}
									}
								}
							}
						}

						newVetos--;
					}
					//

					// Create event-specific histogram and fill it with grouped energies (entries in groupVetoList)
					std::stringstream stringStream;
					stringStream << "event" << numEvents;
					G4String eventName = stringStream.str();
					TH3D* tempH3 = new TH3D(eventName,eventName,iHigh-iLow+1,(iLow-0.9*m/crystalSize)*crystalSize/cm,(iHigh-0.9*m/crystalSize+1)*crystalSize/cm,jHigh-jLow+1,(jLow-0.9*m/crystalSize)*crystalSize/cm,(jHigh-0.9*m/crystalSize+1)*crystalSize/cm,tHigh-tLow+1,tLow*timeRes/ns,(tHigh+1)*timeRes/ns);

					eventList[numEvents-1]->SetLeftXPos(iLow*crystalSize-0.9*m);
					eventList[numEvents-1]->SetBottomYPos(jLow*crystalSize-0.9*m);
					eventList[numEvents-1]->SetStartTime(tLow*timeRes);

					for(int l = 0;l < groupVetoCount;l++)
					{
						double iPos = (groupVetoList[l][0]-0.9*m/crystalSize+0.5)*crystalSize/cm;
						double jPos = (groupVetoList[l][1]-0.9*m/crystalSize+0.5)*crystalSize/cm;
						double time = (groupVetoList[l][2]+.5)*timeRes/ns;
						double content = h3->GetBinContent(groupVetoList[l][0]+1,groupVetoList[l][1]+1,groupVetoList[l][2]+1);
						h3->SetBinContent(groupVetoList[l][0]+1,groupVetoList[l][1]+1,groupVetoList[l][2]+1,0);
						tempH3->Fill(iPos,jPos,time,content);
						eventList[numEvents-1]->AddEnergy(content*MeV);
						delete groupVetoList[l];
					}

					delete groupVetoList;

					eventList[numEvents-1]->SetH3(tempH3);
				}
			}
		}
	}

	for(int l = 0;l < vetoCount;l++)
	{
		delete vetoList[l];
	}

	delete vetoList;

	return numEvents;
}

int** KSADAnalyze::AddToList(int** list,int& count,int i,int j,int t)
{
	int** listOld = list;
	list = new int*[count+1];

	for(int l = 0;l < count;l++)
	{
		list[l] = listOld[l];
	}

	delete listOld;

	list[count] = new int[3];

	list[count][0] = i; list[count][1] = j; list[count][2] = t;

	count++;
	return list;
}

bool KSADAnalyze::InList(int** list,int count,int i,int j,int t)
{
	for(int checkNum = 0;checkNum < count;checkNum++)
	{
		if(list[checkNum][0] == i && list[checkNum][1] == j && list[checkNum][2] == t)
		{
			return true;
		}
	}

	return false;
}

TH3D* KSADAnalyze::GetH3(){return h3;}

void KSADAnalyze::Visualize()
{
	if(h3)
	{
		KSADAnalyzeVis* vis = new KSADAnalyzeVis(this);
		delete vis;
	}
	else
	{
		G4cout << "Please run /analyze/analyze first.\n";
	}
}

KSADAnalyzeEvent** KSADAnalyze::GetEventList()
{
	return eventList;
}

int KSADAnalyze::GetNumEvents()
{
	return numberOfEvents;
}

double KSADAnalyze::GetMaxValue()
{
	return maxValue;
}
